function G = main_func_SE(para)

try
    
    gamma   = 10;
    psi     = 1.5;
    
    main_para
    
    nE    = 501;
    k0    = linspace(k_min,k_max,nK)';
    
    options = optimoptions('fsolve','Display','off');
    beta    = fsolve(@solve_beta,0.99,options,gamma,psi,mu_C,sig_C,P);
    
    main_SDF
    
    lamS    = ones(nS,nK);
    lamD    = ones(nS,nK);
    k       = k0;
    k_bar   = ones(nS,1)*5;
    k_iss   = ones(nS,1)*10;
    
    [kgrid,mgrid]   = meshgrid(k,(1:nS)');
    kgrid           = kgrid(:);
    mgrid           = mgrid(:);
    
    Q_big   = repmat(Q,nK,1); % precompute for speed
    Rf1_mat = repmat(Rf(:,1),1,nK); % [nS,nK]
    
    tic
    main_VFI
    toc

    %if isnan(J), return, end
    
    xvec = [.10 .05 .025 .0125]; % search interval: percent around each target
    pvec = [2,10,50,250]; % increase number of gridpoints in search interval by this factor (relative to the initial, not the previous grid)
    old = [k_def, k_bar, k_iss];
    
    for i = 1:4
        disp(i)
        tic
        Fit_S = griddedInterpolant({1:nS,k},lamS,'linear');
        Fit_D = griddedInterpolant({1:nS,k},lamD,'linear');
        
        x   = xvec(i);
        %         nkf = 200;
        %         kd1 = linspace(k_def(1)*(1-x),k_def(1)*(1+x),nkf)';
        %         kd2 = linspace(k_def(2)*(1-x),k_def(2)*(1+x),nkf)';
        %         kd3 = linspace(k_def(3)*(1-x),k_def(3)*(1+x),nkf)';
        %         kb1 = linspace(k_bar(1)*(1-x),k_bar(1)*(1+x),nkf)';
        %         kb2 = linspace(k_bar(2)*(1-x),k_bar(2)*(1+x),nkf)';
        %         kb3 = linspace(k_bar(3)*(1-x),k_bar(3)*(1+x),nkf)';
        %         ki1 = linspace(k_iss(1)*(1-x),k_iss(1)*(1+x),nkf)';
        %         ki2 = linspace(k_iss(2)*(1-x),k_iss(2)*(1+x),nkf)';
        %         ki3 = linspace(k_iss(3)*(1-x),k_iss(3)*(1+x),nkf)';
        
        
        k_precise = linspace(min(k0),max(k0),length(k0)*pvec(i)-1)'; % subtracting 1 ensures that the original points remain part of the new grid
        k_precise = [k_precise;k_precise(end)+(k_precise(2)-k_precise(1))*(1:10000)']; % add more equally-spaced points above grid (in case target hits kmax)
        keep = false(size(k_precise)); % dummy: all points that fall into one of the search intervals
        keep(k_precise>=(k_def(1)*(1-x)) & k_precise<=(k_def(1)*(1+x))) = true;
        keep(k_precise>=(k_def(2)*(1-x)) & k_precise<=(k_def(2)*(1+x))) = true;
        keep(k_precise>=(k_def(3)*(1-x)) & k_precise<=(k_def(3)*(1+x))) = true;
        keep(k_precise>=(k_def(4)*(1-x)) & k_precise<=(k_def(4)*(1+x))) = true;
        
        keep(k_precise>=(k_bar(1)*(1-x)) & k_precise<=(k_bar(1)*(1+x))) = true;
        keep(k_precise>=(k_bar(2)*(1-x)) & k_precise<=(k_bar(2)*(1+x))) = true;
        keep(k_precise>=(k_bar(3)*(1-x)) & k_precise<=(k_bar(3)*(1+x))) = true;
        keep(k_precise>=(k_bar(4)*(1-x)) & k_precise<=(k_bar(4)*(1+x))) = true;
        
        keep(k_precise>=(k_iss(1)*(1-x)) & k_precise<=(k_iss(1)*(1+x))) = true;
        keep(k_precise>=(k_iss(2)*(1-x)) & k_precise<=(k_iss(2)*(1+x))) = true;
        keep(k_precise>=(k_iss(3)*(1-x)) & k_precise<=(k_iss(3)*(1+x))) = true;
        keep(k_precise>=(k_iss(4)*(1-x)) & k_precise<=(k_iss(4)*(1+x))) = true;
        
        k = [k0;k_precise(keep)]; % original grid + more points in search intervals
        
        
        %k   = [k0;kd1;kd2;kd3;kb1;kb2;kb3;ki1;ki2;ki3];
        k   = sort(unique(k));
        nK  = length(k);
        k_max = max(k);
        
        lamS = Fit_S({1:nS,k});
        lamD = Fit_D({1:nS,k});
        
        [kgrid,mgrid]   = meshgrid(k,(1:nS)');
        kgrid           = kgrid(:);
        mgrid           = mgrid(:);
        
        Q_big   = repmat(Q,nK,1); % precompute for speed
        Rf1_mat = repmat(Rf(:,1),1,nK); % [nS,nK]
        
        main_VFI
        
        new = [k_def, k_bar, k_iss];
        check = abs(new./old-1);
        old = new;
        
        if i==3
            if max(check(:))>xvec(i)
                disp('Warning: Optimal policy outside grid')
            end
        end
        toc
    end
    
    
    main_CDS
    
    money       = ones(nS,nK); % moneyness levels
    main_options
    IV          = NaN(nS,nK,2);
    IV(:,:,1)   = BSIV;
    money       = exp(-BSIV/sqrt(12)); % moneyness levels
    main_options
    IV(:,:,2)   = BSIV;
    IV_ATM      = IV(:,:,1)';
    IV_ATM      = permute(IV_ATM,[2,1]);
    IV_Skew     = IV(:,:,2)'-IV(:,:,1)';
    IV_Skew     = permute(IV_Skew,[2,1]);
    
    main_Ret
    
    main_simu
    
    mom     = [1:2 4:6 9:10 12:15];
    G       = G_m(mom);
    
catch
    
    G   = NaN(1,11);
    
end

